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Abstract We study the quenched dynamics of the momentum distribution of a unitary Bose gas 
under isotropic harmonic confinement within a time-dependent density functional approach based on 
our recently calculated Monte Carlo (MC) bulk equation of state. In our calculations the inter-atomic 
s-wave scattering length of the trapped bosons is suddenly increased to a very large value and the 
real-time evolution of the system is studied. Prompted by the very recent experimental data of 85 Rb 
atoms at unitarity [Nature Phys. 10, 116 (2014)] we focus on the momentum distribution as a function 
of time. Our results suggest that at low momenta, a quasi-stationary momentum distribution is reached 
after a long transient, contrary to what found experimentally for large momenta which equilibrate on 
a time scale shorter than the one for three body losses. 

Keywords Unitary Bose Gas • Time Dependent Density Functional • Local Density approximation • 
Momentum Distribution 


1 Introduction 

Recent experiments claim to have achieved a metastable degenerate gas of ultracold and dilute bosonic 
atoms with infinite s-wave scattering length [l; §;|3; 4}. This is the so-called unitary Bose gas, which is 
characterized by remarkably simple universal laws arising from scale invariance 5]. While the unitary 
Fermi gas has been largely investigated both experimentally and theoretically |6], its bosonic coun¬ 
terpart has been only marginally theoretically addressed [1 ! ! D3 ED; El E| because generally 
considered as experimentally inaccessible M due to the dominant three-body losses at very large 
values of the scattering length. 

By using a Monte Carlo (MC) approach [l3] we have recently investigated the zero-temperature 
properties of a dilute homogeneous Bose gas by tuning the interaction strength of the two-body poten¬ 
tial to achieve arbitrary positive values of the s-wave scattering length a s , while avoiding the formation 
of the self-bound clusters present in the ground-state El • More recently, the MC equation of state has 
been the key ingredient of a local density approximation (LDA) ESI to the energy density functional 
to be used in density functional theory (DFT). Remarkably, the density profiles of a unitary Bose 
gas in a harmonic trap calculated with DFT using such functional compare very well with the MC 
ones. In addition, by using a time-dependent formulation, we have also investigated the excitation 
frequencies as a function of the scattering length. Interestingly, the calculated values for the monopole 
breathing mode, which reproduce the expected limiting values of the ideal and unitary regimes, exhibit 
a non-monotonous behavior as a s varies fl5j . 
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Unfortunately, the experimental data on the unitary Bose gas 0; S i i to compare with are 
quite scarce, and they do not provide an easy tool for validating theoretical approaches. The most 
suitable data in this sense are the time-resolved measurements of the momentum distribution n(k, t) 
of a Bose-condensed gas quenched at unitarity (sudden increase of a s to a very large value) provided 
by the very recent experiment of Makotyn et al. 0|. The main result coming from Ref. 4j is that, 
for k > ks = (67 r 2 n) 1//3 where n is the atom number density, n(k, t) evolves to a quasi-steady-state 
distribution on a time scale shorter than the one characterizing three-body losses. We thus try to 
reproduce the observed phenomenon by solving a nonlinear Schrodinger equation (NLSE) obtained 
from time-dependent DFT (TDDFT) with our LDA functional, based on the MC equation of state 
and where we have introduced a dissipative term to take into account three-body losses, whose role is 
relevant during the quenched dynamics fl 6 j . We must notice however, that we expect our single orbital 
TDDFT to be fully reliable only at low momenta, where the collective long-wavelength dynamics 
dominates. Therefore our analysis will be confined to low momenta, in a range not directly comparable 
with the one studied experimentally. The results we find, are therefore complementary to those of 
Makotyn et al. [|| and suggest that, as expected, while the system equilibrates locally in a short time 
as found by experiments, it takes a very long time to reach quasi equilibrium all over the trap. 


2 Method 


On the basis of the density functional theory D3, a reliable energy functional of the local density n{ r) 
for an inhomogeneous system of interacting bosons at T = 0 is given by 

E[n( r)] = /{*- yj n(r)^ + n(r) s(n(r)) + n(r) t/(r)j d 3 r, (1) 

where the quantum-pressure gradient term takes into account effects due to density variations [18], 
s(n(r)) is the energy per atom of the homogeneous system with density equal to the local density 
and U{ r) describes the external confinement, which we assume to be an isotropic harmonic potential 
U(r) = \mw 2 H r 2 . The values of e(n) have been recently fitted to the results of a MC calculation 03{] for 
a wide range of (positive) values of the scattering length a s characterizing the interparticle interaction. 
In the weakly interacting regime (x = a s /ro -C 1, where r 0 = (3/ (47m)) 1 / 3 is the average distance 
between bosons) the MC results for e(n) are very close to £LHy(n), the universal Bogoliubov prediction 
[Till e b (n) = 5^;(67r 2 n) 2 / 3 as corrected by Lee, Huang and Yang (LHY) [2Cj. In the strong-coupling 
regime (x 1), instead, MC data reach a plateau and, in the unitarity limit (a s —> 00 ), a finite and 
positive energy per particle is found, E/N = 0.70 eg(ra). The equation of state of the homogeneous 
system [Till] from such MC calculation can be well interpolated as: 

/lhy(^) + CI 3 X 3 for x < 0.3 

Cjx 7 + CqX 6 + C 5 X 5 + C 4 X 4 + C 3 X 3 + C 2 X 2 + C\X + cq for 0.3 < x < 0.5 (2) 

b 0 + b\ tanh (b^/x — 1) for x > 0.5 


e B {n) 


with 03 = 0.21, b Q = 0.45, b\ = —0.33, &2 = 0.54, co = 4.75 , ci = —99.72, C 2 = 890.68, C 3 = 
-4309.56, c 4 = 12268.41,c 5 = -20488.00, c 6 = 18568.27 and c 7 = -7052.20 
\ 1/3 


2l|. In ©, /lhyOt) = 


(5&r *[i + 


128 


—£ 3 / 2 ] is the LHY correction to the Bogoliubov prediction. Notice that in the 
deep weak-coupling regime Eq. CD reduces to the familiar Gross-Pitaevskii density functional [22;; 23] 
since e{n) = frjpE^, a s) = 27r?i 2 a s n 2 /m. 

In Ref. [15[ we have shown that the energy functional CD is very accurate in reproducing the 
MC static density profiles of the inhomogeneous unitary Bose gas under harmonic confinement. The 
dynamics of the system can be obtained by generalizing the energy functional m into a density-phase 
action functional A[n(r,t), 6 (r,t)] : where 9(r,t) is the phase of a single-valued quantum-mechanical 
wave function representing the macroscopic wave function of the Bose-Einstein condensate of the 
superfluid [24|, which is related to the superfluid velocity v by the relation v = ( h/m)V9 [25|. Such a 
density-phase action functional is written as 


A[n(r,t),9{r,t)] 


{r[n(r,t),0(r,t)] -£[n(r,f)]j 


dt , 


(3) 
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Fig. 1 Time evolution of the total number N of trapped 85 Rb atoms in the unitary Bose gas. Solid line: 
numerical integration of the dissipative nonlinear Schrodinger equation 0. Filled squares: experimental data 
from Ref. [4|. 


where 

T[n(r,t),9(r,t)} = j j-n(r,t) ^ ( V6> M)) 2 ) } d 3 * ( 4 ) 

is the kinetic Lagrangian of Popov |26], and E[n(r,t)] is given by Eq. © under the assumption of a 
time-dependent local density. It is straightforward to derive the equations of superfluid hydrodynamics 
with a quantum-pressure term (directly related to the von Weizsacker gradient term of Eq. (0)) by 
extremizing the action functional (0) |27j |. Moreover, introducing the wave function 

tf(r,i) = e ie ^ (5) 


the equations of superfluid hydrodynamics can be re-written in terms of a nonlinear Schrodinger 
equation 


ih 


d\P(r, t) 
dt 


' h 2 v 2 

2 TO 


+ U (r) + 


d{ne) 

dn 


^(M) , 


( 6 ) 


which is the Euler-Lagrange equation obtained from the action functional (0) taking into account Eq. 

©• 

We stress that Eq. © can be alternatively obtained from time dependent density functional theory 
(TDDFT) [13 123121| in the Local Density Approximation using a single Kohn-Sham orbital to describe 
the degenerate boson system [13;Ell • The computational approach based on E q. © has been adopted, 
for instance, to describe superfluid helium j30t f3l'lj , ultracold bosonic atoms (27; H. and also superfluid 
fermions in the BCS-BEC crossover (32: l33l|7 

In the dynamics of the unitary Bose gas three-body losses play a relevant dissipative role, especially 
close or at unitarity. As done in previous applications of Eq. © [ll|, we model the effect of three-body 
losses by adding a phenomenological dissipative term — (with L 3 = 9 ■ 10 -23 cm 6 /sec [i|) in 
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Fig. 2 Scaled average radius (r 2 ) 1 / 2 / an of the 85 Rb atomic cloud as a function of time t. Solid line: numerical 
integration of the dissipative NLSE, Eq. 0. Dotted line: numerical integration of the dissipative GPE, i.e. 
Eq. 0 with e(n) = £gpe (n). In both cases the scattering length a B is very large but finite: a s = 5 • 10 5 a o, 
with a o = 0.53 • 10~ 10 m the Bohr radius. The two horizontal dashed lines give the experimentally-estimated 
average radius of the cloud up to the time t = 500 fis [J]. an = \fh/(muiH) is the characteristic length of the 
harmonic confinement. 


equation (|5jl. Thus, the dissipative NLSE we use for our numerical simulations is given by 

ihL 3 


,n — t ~ 


" 2v! + U( r) + 8( ” e) 


2 m 


dn 


-l^r ,t)r 


^(r ,t) . 


(7) 


We have numerically solved Eq. 0 to obtain the real-time evolution closely simulating the ex¬ 
perimental conditions of Makotyn et al. Q. Therefore we consider as initial configuration a cloud 
of N = 70 0 0 0 85 i ?6 atoms, confined in a spherical harmonic trap with frequency loh = 10 Hz and 
scattering length a s = 150 ao (where ao is the Bohr radius) prepared in its ground state by evolving 
in imaginary time Eq. 0 without the dissipative term (i.e. with L 3 = 0). Then we switch on the 
dissipative term (with L 3 = 9 ■ 10~ 23 cm 6 /sec), increase the scattering length to a very high, but 
otherwise arbitrary, value a s = 5 x 10 5 a o, and let the system evolve in real time for a time t. At such 
time we analyze the momentum distribution of the particles, as described in the following. Notice that 
according to our calculated MC equation of state 0, the chosen value for a s is practically equivalent 
to a s = Too, i.e. it gives the same energy per atom as in the truly unitary limit [13|. 


3 Numerical results 

Due to the presence of the three-body recombination, whose rate increases as of, atoms are continuously 
disappearing from the trap. At early times, t < 350 fi s, experimental data are compatible with an 
exponential decay with a time constant of about 630 /xs [4j . In our simulations the number of particles 
N in the trap suffers a similar depletion, driven by the phenomenological dissipative term proportional 
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Fig. 3 Momentum distributions n(k ) of the bosonic cloud of 85 Rb atoms obtained at increasing time steps by 
solving Eq. ©■ The narrowest distribution corresponds to t = 0 while the broadest one to t = 1000 /rs. There 
are 37 distributions obtained at increasing time steps At = 27.02 /.ts. The data are shown as a function of the 
reduced wave-vector k = k/ks, where fcs = ( 67 r 2 n ) 1//3 and n is the density at the center of the trap. 


to L 3 . The time-dependent behavior of the total number of trapped atoms can be monitored by solving 
Eq- (0 and by computing 

N(t) = j |^(r,f )| 2 d 3 r . ( 8 ) 

The obtained N(t) is shown as a solid line in Fig.|T]and is compared to the experimental data of Ref. Q, 
reported as filled squares. Our DFT results are fairly comparable with the experimental points, where 
these are available, but we note that for later times our decay is far from being exponential. It seems to 
be characterized by two different time scales: one, fast, driving the first depletion of the atoms within 
the trap, that covers the experimental data range; and a second one, which is much slower, dictating 
the emptying of the trap at very large times. We find that, in this second regime, a relevant number 
of atoms still populate the trap, and the depletion is much slower, leading to an apparent stationary 
condition, as discussed in the following, and allowing for measurements of properties that give the 
impression of being equilibrated. 

The size of trapped cloud predicted on the basis of our DFT calculation is also in reasonable 
agreement with the experimental data, given their uncertainty. In Fig.[2]we report the average radius of 
the bosonic cloud as a function of time t. The solid line refers to the results obtained using the dissipative 
NLSE 0 with the MC equation of state © and as = 5 x 10 5 ao- They are fully compatible with the 
experimental findings pj reported as horizontal dashed lines in Fig. [3 For the sake of comparison, 
in Fig. [2] we also plot with a dotted line the results for the same average radius, computed using 
the dissipative NLSE Eq. [G] with the same finite scattering length a s and the same three-body loss 
coefficient L 3 as above, but with the energy density e(n) from the Gross-Pitaevskii equation of state 
instead than from MC. The figure clearly shows that the radius obtained within the GPE theory 
(dotted line) rapidly exits the experimental range as measured in the first 500/iS Q and displays local 
oscillations due to interference effects in the outer part of the atomic cloud interacting with the trap 
walls. 
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t (/i-s) 


Fig. 4 Time evolution of the number of atoms N ( ki ) with momentum ki (expressed in 
The 5 curves correspond to 5 values of ki : 0.08 (solid line), 0.12 (dotted line), 0.16 
dashed line), 0.24 (dot-dashed line). The curves have been arbitrarily normalized to 
0.08) (solid line). 


units of ks = 

(dashed line), 0.20 (long 
the maximum of N(ki = 


Given the solution ^(r, i) of the dissipative NLSE 0, the time dependent momentum distribution 
is easily calculated as: 


i(k,f) = 


J'\P(r,t ) exp(ik-r) d 3 r 


(9) 


In Fig.[3]we plot n(k, t) obtained at different times, up to a maximum value of 1 ms, with k = k/ks , 
where ks = ( 67 r 2 ?!) 1 / 3 and n is the density at the center of the trap. Notice that in the figure we report 
n(k,t) at low momenta (0 < k < 0.3), because our single-orbital TDDFT is supposed to be fully 
reliable only at long wavelengths. During the first 300 /zs we observe an expansion of the cloud both in 
real and in momentum spaces, and then a quasi-stationary momentum distribution seems to develop 
for larger times. 

However, we must point out that at longer times (not shown in the figure) n(k,t) is still changing. 
This could be ascribed to the fact that, due to the sudden quench of a s to so high values, the bosonic 
cloud expands until the wave-front interacts with the steeper part of the trap walls and is reflected 
back, letting the cloud to shrink again at much larger times. Thus the momentum distribution first 
decreases (close to the turning point) and then it increases again as the contraction of the cloud occurs. 

To better clarify the time scales involved in the momentum changes, we plot in Fig. [3] the time 
evolution of the number of atoms N{k{) with a given momentum ki for different values of ki during 
the first ~ 1 ms after the sudden quench. We note that the evolution of the distribution for these small 
momenta is not uniform, being oscillatory first, with the amplitude of the oscillations extinguishing on 
a time scale of 700 — 800 p, s. 

This is to contrasted with what found in the experiments of Ref.|4] for higher (k > ks) components 
in the momentum distribution. There, a shorter time scale governs the saturation of the k-component, 
increasing from about 100/xs to about 300 — 400 fis as k decreases from k/ks ~ 1.2 —1.3 to k/ks ~ 0.8. 
Our results complement these observations, in finding a still larger time-scale when k/ks ~ 0.1 — 0.3. 
Our results can be compared also with recent studies on the dynamics of n{k , t) after a sudden quench 
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to very large value of a s 0. In particular, calculations based on a bath approach (34j lead to a 
multistep equilibration process, where larger k modes equilibrate faster (having a shorter relaxation 
time), the oscillations in n(k) are damped on intermediate times, leading to an apparent equilibrated 
state, and the full equilibration is reached only on a much larger time scale, when also the condensate 
attains its final time-independent state. 


4 Conclusions 

We have numerically investigated the quenched dynamics of a Bose gas of 85 Rb atoms under isotropic 
harmonic confinement after the sudden increase of the s-wave scattering length from a s = 150 ao to 

5 • 10 5 a o- To take into account three body losses we have introduced a dissipative term into the time 
dependent nonlinear Schrodinger equation (TDNLSE) obtained from a local approximation of a TD 
energy functional. We have shown that, while the equation obtained by using the Gross-Pitaevskii 
equation of state in the TD energy functional is unable to fully capture the experimental quenched 
dynamics of the atomic cloud, the corresponding equation derived from a functional containing the 
energy density fitted to our recently calculated Monte Carlo bulk equation of state [l3| gives values of 
the average radius of the cloud in fairly good agreement with very recent experimental data Q. We 
take this result as supporting the use of our MC equation of state in the energy density functional. 

The solution of our dissipative TDNLSE shows a fast depletion of the condensate at short times, 
t < 300 [is due to three-body losses, accompanied by the relative increase in the population of large 
momenta. Also this feature is in agreement with the experimental findings Q. As shown in Fig. Q] 
though, the momentum densities at various momenta reach a ’quasi equilibrium’ distribution, at very 
large times of about 1000 /is. One must notice however that even the larger momenta reported in 
our figures are much smaller than those studied experimentally, corresponding to values for which our 
TDDFT is questionable. We reconcile our results with those reported in Ref. [|J by observing that 
they refer to different length scales: while from the experiment a fast equilibration on a short length 
scale has been found, reminiscent of local equilibrium, our calculation shows that a much longer time 
is needed to get quasi-equilibrium on a large length scale. The authors of Ref. @] interpret their results 
(see their fig. 4) as a saturation of the value of n(k ) at times of the order of 200 p,s, shorter than the 
time scale for three body losses (300 /zs), thus indicating a different mechanism as responsible for the 
equilibration of the momentum distribution. Such a mechanism could be the interaction between the 
condensed and non condensed fractions of boson atoms which is not present in our equation, since the 
dissipative term in our equation lets the total number of atoms to decrease while leaving the system 
in a degenerate state. 

The effect of the interaction between the condensate and excited atoms on the dynamics of mo¬ 
mentum distribution of a Bose gas after a rapid quench to a very large value of a s has been explicitly 
addressed at in some recent papers [34]; [33 [33 [33; [3g |39| and supports the existence of different times 
domains in the dynamics of momentum distribution. In Ref. 1341 . in particular, they correspond to a 
fast depletion of the condensate followed or accompanied by a slow re-adjustment of the momentum 
distribution and, only at very large times, the reaching of a metastable state. We will explore the 
possibility of coupling the wave function obeying our TDNLSE to a bath [34] describing the small 
noncondensed fraction with the goal of better reproducing the experimental data at large momenta 
and short times. 

The authors are grateful to J. P. Corson for a critical reading of a previous version of this report and 
for useful remarks. The authors acknowledge for partial support the Ministero Istruzione Universita 
Ricerca (PRIN project 2010LLKJBX). 
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